using namespace std;
#include<iostream>
#include<iomanip>
#include<fstream>
#include<vector>
#include<stdio.h>
#include<stdlib.h>
#include<map>
#include<set>
#include"BA2.h"
#include<math.h>
#include<gsl/gsl_rng.h>
#include<gsl/gsl_randist.h>
#include<gsl/gsl_sf_gamma.h>
#include"cpptime.h"

extern std::vector<string> cmdList;
extern std::vector<Indiv> sample;
extern std::map<string,int> popToindex;
extern std::map<int,string> indexTopop;
extern std::map<string,int> labelToindex;
extern std::map<int,string> indexTolabel;
extern std::set<string> labelSet;
extern std::set<string> popSet;
extern std::vector<long int> popCounts;
extern std::vector< vector< vector<long int> > > alleleCounts;
extern std::vector< set<string> > alleleSet;
extern std::vector< map<string,int> > alleleMap;
extern std::vector< map<int,string> > intToallele;
extern std::set<string> locusSet;
extern std::map<int,string> indexTolocus;
extern std::map<string,int> locusToindex;
extern double maxMigration;
extern bool infileOpen;

// run MCMC without data
// #define NODATA

// run MCMC without modifying allele freqs
// #define NOALLELEFREQ

// run MCMC without modifying migration rates
// #define NOMIGRATIONRATE

// run MCMC without migrant count prior
// #define NOMIGRANTCOUNTPRIOR

// run MCMC without modifying inbreeding coefficient
// #define NOINBREEDING

// use uniform allele freqs to initialize MCMC (otherwise use empirical freqs)
// #define UNIFORM

// constructor for BayesAss class
BayesAss::BayesAss()
{
  for(int j=0;j<static_cast<int>(sample.size());j++)
    {
      AugmentIndiv tempIndiv;
      tempIndiv.migrantPopulation=popToindex[sample[j].sampledPopulation];
      tempIndiv.ancestryAge=0;
      tempIndiv.fstat=0.0;
      augmentedSample.push_back(tempIndiv);
    }
  vector< vector<double> > tempPopln;
  for(int j=0;j<=static_cast<int>(sample[0].genotypes.size())-1;j++)
    {
      vector<double> tempLocus;
      for(int k=0; k<=static_cast<int>(alleleSet[j].size())-1; k++)
	tempLocus.push_back(0);
      tempPopln.push_back(tempLocus);
    }
  // create structure for allele freqs
  for(int i=0;i<=static_cast<int>(popSet.size())-1;i++)
    {
      alleleFreqs.push_back(tempPopln);
    }

#ifdef UNIFORM
  // initialize allele freqs to uniform
  for(int i=0;i<static_cast<int>(alleleFreqs.size());i++)
    for(int j=0;static_cast<int>(j<static_cast<int>(alleleFreqs[i].size()));j++)
      {
	// set freq for missing data to 1 (e.g., probability=1)
	alleleFreqs[i][j][0]=1.0;
	// set remaining allele freqs to be uniform 
	for(int k=1;k<static_cast<int>(alleleFreqs[i][j].size());k++)
	  alleleFreqs[i][j][k]=1.0/static_cast<int>(alleleFreqs[i][j].size()-1);
      }
#endif

#ifndef UNIFORM
  // initialize allele freqs to empirical freqs
  for(int i=0;i<static_cast<int>(alleleFreqs.size());i++)
    for(int j=0;static_cast<int>(j<static_cast<int>(alleleFreqs[i].size()));j++)
      {
	// set freq for missing data to 1 (e.g., probability=1)
	alleleFreqs[i][j][0]=1.0;
	double sumAFs=0.0;
	for(int k=1;k<static_cast<int>(alleleFreqs[i][j].size());k++)
	  sumAFs+=alleleCounts[i][j][k];
	for(int k=1;k<static_cast<int>(alleleFreqs[i][j].size());k++)
	  alleleFreqs[i][j][k]=(alleleCounts[i][j][k]+0.00000001)/sumAFs;
      }
#endif

  // create structure for average posterior allele freqs
  avgAlleleFreqs=alleleFreqs;

  // create structure for migration rates and average posterior migration rates
  for(int i=0;i<static_cast<int>(popSet.size());i++)
    {
      vector<double> tempRow;
      for(int j=0;j<static_cast<int>(popSet.size());j++)
	{
	  if(i==j) 
	    tempRow.push_back(0.95);
	  else
	    tempRow.push_back(0.05/(static_cast<int>(popSet.size())-1.0));
	}
      migrationRates.push_back(tempRow);
    }
  avgMigrationRates=migrationRates;
  m2MigrationRates=migrationRates;

  // create structure for migrant counts
  vector<long int> tempRow1;
  vector<vector<long int> > tempRow2; 
  for(int k=0;k<=2;k++)
    tempRow1.push_back(0);
  for(int i=0;i<static_cast<int>(popSet.size());i++)
    tempRow2.push_back(tempRow1);
  for(int i=0;i<static_cast<int>(popSet.size());i++)
	migrantCounts.push_back(tempRow2);

  // intialize migrant counts
  for(int i=0;i<static_cast<int>(sample.size());i++)
    ++migrantCounts[popToindex[sample[i].sampledPopulation]][augmentedSample[i].migrantPopulation][augmentedSample[i].ancestryAge];
}

double BayesAss::slidingWindowProposal(double perturb, double lowerBound, double upperBound, double origValue)
{
double x=origValue+perturb;
 while((x>upperBound)||(x<lowerBound))
   {
     if(x>upperBound) x=upperBound-fabs(upperBound-x);
     else x=lowerBound+fabs(lowerBound-x);
   }
  return(x);
}

// functions for calculating probability of migrant counts given migration rates

// global probability
double BayesAss::BAmigratelogP(Parameters pvalues)
{
  double logPr=0.0;
  vector<double> sumOffDiagonalRates;
  for(int i=0;i<static_cast<int>(popSet.size());i++)  
    {
      double tempSum=0.0;
      for(int j=0;j<static_cast<int>(popSet.size());j++)
	{
	  if(i!=j) tempSum+=migrationRates[j][i];
	}
      sumOffDiagonalRates.push_back(tempSum);
    }

  double myResult=0.0;
  
  // if Ngens==1 use 1st and 2nd generation migrant ancestries
  if(pvalues.Ngens==1)
    {
      for(int i=0;i<static_cast<int>(popSet.size());i++)
	for(int j=0;j<static_cast<int>(popSet.size());j++)
	  {
	    for(int k=0;k<=2;k++)
	      {
		if((k==0)&&(i==j))
		  logPr+=(migrantCounts[i][j][k]*log(1.0-3.0*sumOffDiagonalRates[i])-
			  gsl_sf_lngamma(migrantCounts[i][j][k]+1.0));
		else 
		  if((k==1)&&(i!=j))
		    logPr+=(migrantCounts[i][j][k]*log(migrationRates[j][i])-
			    gsl_sf_lngamma(migrantCounts[i][j][k]+1.0));
		  else
		    if((k==2)&&(i!=j))
		      logPr+=(migrantCounts[i][j][k]*log(2.0*migrationRates[j][i])-
			      gsl_sf_lngamma(migrantCounts[i][j][k]+1.0));	      
	      }	
	  }
      for(int i=0;i<static_cast<int>(popSet.size());i++)
	myResult+=gsl_sf_lngamma(popCounts[i]+1.0);
    }
  // else if Ngens=0 use only 1st generation migrant ancestries
  else
    {
      for(int i=0;i<static_cast<int>(popSet.size());i++)
	for(int j=0;j<static_cast<int>(popSet.size());j++)
	  {
	    for(int k=0;k<=2;k++)
	      {
		if((k==0)||(k==2))
		  logPr+=(migrantCounts[i][j][k]*log(1.0-sumOffDiagonalRates[i])-
			  gsl_sf_lngamma(migrantCounts[i][j][k]+1.0));
		else 
		  if((k==1)&&(i!=j))
		    logPr+=(migrantCounts[i][j][k]*log(migrationRates[j][i])-
			    gsl_sf_lngamma(migrantCounts[i][j][k]+1.0));
	      }
	  }	
    }
  return(logPr+myResult);
}

// updated probability when new migration rates are proposed
double BayesAss::BAmigratelogP(Parameters pvalues, vector< vector <double> > proposedMigrationRates)
{
  double logPr=0.0;
  vector<double> sumOffDiagonalRates;
  for(int i=0;i<static_cast<int>(popSet.size());i++)  
    {
      double tempSum=0.0;
      for(int j=0;j<static_cast<int>(popSet.size());j++)
	{
	  if(i!=j) tempSum+=proposedMigrationRates[j][i];
	}
      sumOffDiagonalRates.push_back(tempSum);
    }

  double myResult=0.0;
  
  // if Ngens==1 use 1st and 2nd generation migrant ancestries
  if(pvalues.Ngens==1)
    {
      for(int i=0;i<static_cast<int>(popSet.size());i++)
	for(int j=0;j<static_cast<int>(popSet.size());j++)
	  for(int k=0;k<=2;k++)
	    {
	      if(k==0)
		logPr+=(migrantCounts[i][j][k]*log(1.0-3.0*sumOffDiagonalRates[i])-
			gsl_sf_lngamma(migrantCounts[i][j][k]+1.0));
	      else 
		if(k==1)
		  logPr+=(migrantCounts[i][j][k]*log(proposedMigrationRates[j][i])-
			  gsl_sf_lngamma(migrantCounts[i][j][k]+1.0));
		else
		  if(k==2)
		    logPr+=(migrantCounts[i][j][k]*log(2.0*proposedMigrationRates[j][i])-
			    gsl_sf_lngamma(migrantCounts[i][j][k]+1.0));	      
	    }
      for(int i=0;i<static_cast<int>(popSet.size());i++)
	myResult+=gsl_sf_lngamma(popCounts[i]+1.0);
    }
  // else if Ngens=0 use only 1st generation migrant ancestries
  else
    {
      for(int i=0;i<static_cast<int>(popSet.size());i++)
	for(int j=0;j<static_cast<int>(popSet.size());j++)
	  {
	    for(int k=0;k<=2;k++)
	      {
		if((k==0)||(k==2))
		  logPr+=(migrantCounts[i][j][k]*log(1.0-sumOffDiagonalRates[i])-
			  gsl_sf_lngamma(migrantCounts[i][j][k]+1.0));
		else 
		  if((k==1)&&(i!=j))
		    logPr+=(migrantCounts[i][j][k]*log(proposedMigrationRates[j][i])-
			    gsl_sf_lngamma(migrantCounts[i][j][k]+1.0));
	      }
	  }	

    }
  return(logPr+myResult);
}

// updated probability when new migrant counts (individual migrant ancestries) are proposed
double BayesAss::BAmigratelogP(Parameters pvalues, vector< vector< vector<long int> > > proposedMigrantCounts)
{
  double logPr=0.0;
  vector<double> sumOffDiagonalRates;
  for(int i=0;i<static_cast<int>(popSet.size());i++)  
    {
      double tempSum=0.0;
      for(int j=0;j<static_cast<int>(popSet.size());j++)
	{
	  if(i!=j) tempSum+=migrationRates[j][i];
	}
      sumOffDiagonalRates.push_back(tempSum);
    }

  double myResult=0.0;
  
  // if Ngens==1 use 1st and 2nd generation migrant ancestries
  if(pvalues.Ngens==1)
    {
      for(int i=0;i<static_cast<int>(popSet.size());i++)
	for(int j=0;j<static_cast<int>(popSet.size());j++)
	  for(int k=0;k<=2;k++)
	    {
	      if(k==0)
		logPr+=(proposedMigrantCounts[i][j][k]*log(1.0-3.0*sumOffDiagonalRates[i])-
			gsl_sf_lngamma(proposedMigrantCounts[i][j][k]+1.0));
	      else 
		if(k==1)
		  logPr+=(proposedMigrantCounts[i][j][k]*log(migrationRates[j][i])-
			  gsl_sf_lngamma(proposedMigrantCounts[i][j][k]+1.0));
		else
		  if(k==2)
		    logPr+=(proposedMigrantCounts[i][j][k]*log(2.0*migrationRates[j][i])-
			    gsl_sf_lngamma(proposedMigrantCounts[i][j][k]+1.0));	      
	    }
      for(int i=0;i<static_cast<int>(popSet.size());i++)
	myResult+=gsl_sf_lngamma(popCounts[i]+1.0);
    }
  // else if Ngens=0 use only 1st generation migrant ancestries
  else
    {
      for(int i=0;i<static_cast<int>(popSet.size());i++)
	for(int j=0;j<static_cast<int>(popSet.size());j++)
	  {
	    for(int k=0;k<=2;k++)
	      {
		if((k==0)||(k==2))
		  logPr+=(proposedMigrantCounts[i][j][k]*log(1.0-sumOffDiagonalRates[i])-
			  gsl_sf_lngamma(proposedMigrantCounts[i][j][k]+1.0));
		else 
		  if((k==1)&&(i!=j))
		    logPr+=(proposedMigrantCounts[i][j][k]*log(migrationRates[j][i])-
			    gsl_sf_lngamma(proposedMigrantCounts[i][j][k]+1.0));
	      }
	  }	

    }
  return(logPr+myResult);
}

// functions for calculating log likelihood of genotypes given population ancestries

// calculate logL for indiv (all loci)
double BayesAss::BAlocallogL(int indiv) 
{
  double logPr=0.0;
  if((augmentedSample[indiv].ancestryAge==0)||(augmentedSample[indiv].ancestryAge==1))
    { 
      for(int j=0;j<static_cast<int>(sample[indiv].genotypes.size());j++)
	{
	  if(sample[indiv].genotypes[j].g1==sample[indiv].genotypes[j].g2) // if homozygous
	    logPr+=log((1.0-augmentedSample[indiv].fstat)*alleleFreqs[augmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g2]]*
		       alleleFreqs[augmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g2]]+
		       augmentedSample[indiv].fstat*alleleFreqs[augmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g2]]);
	  if(sample[indiv].genotypes[j].g1!=sample[indiv].genotypes[j].g2) // if heterozgous
	    logPr+=(log(2.0)+log(1.0-augmentedSample[indiv].fstat)+
		    log(alleleFreqs[augmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g1]])+
		    log(alleleFreqs[augmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g2]]));
	}
    }
  else
    if(augmentedSample[indiv].ancestryAge==2)
      {
      for(int j=0;j<static_cast<int>(sample[indiv].genotypes.size());j++)
	{
	  if(sample[indiv].genotypes[j].g1==sample[indiv].genotypes[j].g2) // if homozygous
	    logPr+=log(alleleFreqs[augmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g1]])+
	      log(alleleFreqs[popToindex[sample[indiv].sampledPopulation]][j][alleleMap[j][sample[indiv].genotypes[j].g1]]);
	   if(sample[indiv].genotypes[j].g1!=sample[indiv].genotypes[j].g2) // if heterozgous
	    logPr+=log(alleleFreqs[augmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g1]]*
		       alleleFreqs[popToindex[sample[indiv].sampledPopulation]][j][alleleMap[j][sample[indiv].genotypes[j].g2]]+
		       alleleFreqs[augmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g2]]*
		       alleleFreqs[popToindex[sample[indiv].sampledPopulation]][j][alleleMap[j][sample[indiv].genotypes[j].g1]]);
	}
      }
  return(logPr);
}

// calculate logL for indiv (all loci) with proposed change in migrant ancestry
double BayesAss::BAlocallogL(int indiv,  vector<AugmentIndiv> proposedAugmentedSample) 
{
  double logPr=0.0;
  if((proposedAugmentedSample[indiv].ancestryAge==0)||(proposedAugmentedSample[indiv].ancestryAge==1))
    { 
      for(int j=0;j<static_cast<int>(sample[indiv].genotypes.size());j++)
	{
	  if(sample[indiv].genotypes[j].g1==sample[indiv].genotypes[j].g2) // if homozygous
	    logPr+=log((1.0-proposedAugmentedSample[indiv].fstat)*alleleFreqs[proposedAugmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g2]]*
		       alleleFreqs[proposedAugmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g2]]+
		       proposedAugmentedSample[indiv].fstat*alleleFreqs[proposedAugmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g2]]);
	  if(sample[indiv].genotypes[j].g1!=sample[indiv].genotypes[j].g2) // if heterozgous
	    logPr+=(log(2.0)+log(1.0-proposedAugmentedSample[indiv].fstat)+
		    log(alleleFreqs[proposedAugmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g1]])+
		    log(alleleFreqs[proposedAugmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g2]]));
	}
    }
  else
    if(proposedAugmentedSample[indiv].ancestryAge==2)
      {
      for(int j=0;j<static_cast<int>(sample[indiv].genotypes.size());j++)
	{
	  if(sample[indiv].genotypes[j].g1==sample[indiv].genotypes[j].g2) // if homozygous
	    logPr+=log(alleleFreqs[proposedAugmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g1]])+
	      log(alleleFreqs[popToindex[sample[indiv].sampledPopulation]][j][alleleMap[j][sample[indiv].genotypes[j].g1]]);
	   if(sample[indiv].genotypes[j].g1!=sample[indiv].genotypes[j].g2) // if heterozgous
	    logPr+=log(alleleFreqs[proposedAugmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g1]]*
		       alleleFreqs[popToindex[sample[indiv].sampledPopulation]][j][alleleMap[j][sample[indiv].genotypes[j].g2]]+
		       alleleFreqs[proposedAugmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g2]]*
		       alleleFreqs[popToindex[sample[indiv].sampledPopulation]][j][alleleMap[j][sample[indiv].genotypes[j].g1]]);
	}
      }
  return(logPr);
}

// calculate logL for indiv (all loci) with proposed change in allele freq
double BayesAss::BAlocallogL(int indiv, vector< vector< vector<double> > > modAlleleFreqs) 
{
  double logPr=0.0;
  if((augmentedSample[indiv].ancestryAge==0)||(augmentedSample[indiv].ancestryAge==1))
    { 
      for(int j=0;j<static_cast<int>(sample[indiv].genotypes.size());j++)
	{
	  if(sample[indiv].genotypes[j].g1==sample[indiv].genotypes[j].g2) // if homozygous
	    logPr+=log((1.0-augmentedSample[indiv].fstat)*modAlleleFreqs[augmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g2]]*
		       modAlleleFreqs[augmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g2]]+
		       augmentedSample[indiv].fstat*modAlleleFreqs[augmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g2]]);
	  if(sample[indiv].genotypes[j].g1!=sample[indiv].genotypes[j].g2) // if heterozgous
	    logPr+=(log(2.0)+log(1.0-augmentedSample[indiv].fstat)+
		    log(modAlleleFreqs[augmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g1]])+
		    log(modAlleleFreqs[augmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g2]]));
	}
    }
  else
    if(augmentedSample[indiv].ancestryAge==2)
      {
      for(int j=0;j<static_cast<int>(sample[indiv].genotypes.size());j++)
	{
	  if(sample[indiv].genotypes[j].g1==sample[indiv].genotypes[j].g2) // if homozygous
	    logPr+=log(modAlleleFreqs[augmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g1]])+
	      log(modAlleleFreqs[popToindex[sample[indiv].sampledPopulation]][j][alleleMap[j][sample[indiv].genotypes[j].g1]]);
	   if(sample[indiv].genotypes[j].g1!=sample[indiv].genotypes[j].g2) // if heterozgous
	    logPr+=log(modAlleleFreqs[augmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g1]]*
	      modAlleleFreqs[popToindex[sample[indiv].sampledPopulation]][j][alleleMap[j][sample[indiv].genotypes[j].g2]]+
	      modAlleleFreqs[augmentedSample[indiv].migrantPopulation][j][alleleMap[j][sample[indiv].genotypes[j].g2]]*
	      modAlleleFreqs[popToindex[sample[indiv].sampledPopulation]][j][alleleMap[j][sample[indiv].genotypes[j].g1]]);
	}
      }
  return(logPr);
}

// calculate logL for indiv (single locus)
double BayesAss::BAlocallogL(int indiv, int aLocus) 
{
  double logPr=0.0;
  if((augmentedSample[indiv].ancestryAge==0)||(augmentedSample[indiv].ancestryAge==1))
    { 
      if(sample[indiv].genotypes[aLocus].g1==sample[indiv].genotypes[aLocus].g2) // if homozygous
	    logPr+=log((1.0-augmentedSample[indiv].fstat)*alleleFreqs[augmentedSample[indiv].migrantPopulation][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g2]]*
		       alleleFreqs[augmentedSample[indiv].migrantPopulation][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g2]]+
		       augmentedSample[indiv].fstat*alleleFreqs[augmentedSample[indiv].migrantPopulation][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g2]]);
	  if(sample[indiv].genotypes[aLocus].g1!=sample[indiv].genotypes[aLocus].g2) // if heterozgous
	    logPr+=(log(2.0)+log(1.0-augmentedSample[indiv].fstat)+
		    log(alleleFreqs[augmentedSample[indiv].migrantPopulation][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g1]])+
		    log(alleleFreqs[augmentedSample[indiv].migrantPopulation][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g2]]));
    }
  else
    if(augmentedSample[indiv].ancestryAge==2)
      {
	  if(sample[indiv].genotypes[aLocus].g1==sample[indiv].genotypes[aLocus].g2) // if homozygous
	    logPr+=log(alleleFreqs[augmentedSample[indiv].migrantPopulation][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g1]])+
	      log(alleleFreqs[popToindex[sample[indiv].sampledPopulation]][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g1]]);
	   if(sample[indiv].genotypes[aLocus].g1!=sample[indiv].genotypes[aLocus].g2) // if heterozgous
	    logPr+=log(alleleFreqs[augmentedSample[indiv].migrantPopulation][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g1]]*
		       alleleFreqs[popToindex[sample[indiv].sampledPopulation]][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g2]]+
		       alleleFreqs[augmentedSample[indiv].migrantPopulation][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g2]]*
		       alleleFreqs[popToindex[sample[indiv].sampledPopulation]][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g1]]);
	}
  return(logPr);
}

// calculate logL for indiv (single locus) with proposed change in allele freqs
double BayesAss::BAlocallogL(int indiv, int aLocus, vector< vector< vector<double> > > modAlleleFreqs) 
{
  double logPr=0.0;
  if((augmentedSample[indiv].ancestryAge==0)||(augmentedSample[indiv].ancestryAge==1))
    { 
      if(sample[indiv].genotypes[aLocus].g1==sample[indiv].genotypes[aLocus].g2) // if homozygous
	    logPr+=log((1.0-augmentedSample[indiv].fstat)*modAlleleFreqs[augmentedSample[indiv].migrantPopulation][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g2]]*
		       modAlleleFreqs[augmentedSample[indiv].migrantPopulation][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g2]]+
		       augmentedSample[indiv].fstat*modAlleleFreqs[augmentedSample[indiv].migrantPopulation][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g2]]);
	  if(sample[indiv].genotypes[aLocus].g1!=sample[indiv].genotypes[aLocus].g2) // if heterozgous
	    logPr+=(log(2.0)+log(1.0-augmentedSample[indiv].fstat)+
		    log(modAlleleFreqs[augmentedSample[indiv].migrantPopulation][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g1]])+
		    log(modAlleleFreqs[augmentedSample[indiv].migrantPopulation][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g2]]));
    }
  else
    if(augmentedSample[indiv].ancestryAge==2)
      {
	  if(sample[indiv].genotypes[aLocus].g1==sample[indiv].genotypes[aLocus].g2) // if homozygous
	    logPr+=log(modAlleleFreqs[augmentedSample[indiv].migrantPopulation][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g1]])+
	      log(modAlleleFreqs[popToindex[sample[indiv].sampledPopulation]][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g1]]);
	   if(sample[indiv].genotypes[aLocus].g1!=sample[indiv].genotypes[aLocus].g2) // if heterozgous
	    logPr+=log(modAlleleFreqs[augmentedSample[indiv].migrantPopulation][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g1]]*
		       modAlleleFreqs[popToindex[sample[indiv].sampledPopulation]][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g2]]+
		       modAlleleFreqs[augmentedSample[indiv].migrantPopulation][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g2]]*
		       modAlleleFreqs[popToindex[sample[indiv].sampledPopulation]][aLocus][alleleMap[aLocus][sample[indiv].genotypes[aLocus].g1]]);
	}
  return(logPr);
}

// calculate global logL (all indiv, all loci)
double BayesAss::BAgloballogL()
{
  double logPr=0.0;
  for(int i=0;i<static_cast<int>(sample.size());i++)
    logPr+=BAlocallogL(i);
  return(logPr);
}

// calculate global logL (all indiv, single locus)
double BayesAss::BAgloballogL(int aLocus)
{
  double logPr=0.0;
  for(int i=0;i<static_cast<int>(sample.size());i++)
    logPr+=BAlocallogL(i,aLocus);
  return(logPr);
}

// calculate global logL (all indiv, single locus) with proposed change in allele freqs
double BayesAss::BAgloballogL(int aLocus, vector< vector< vector<double> > > modAlleleFreqs)
{
  double logPr=0.0;
  for(int i=0;i<static_cast<int>(sample.size());i++)
    logPr+=BAlocallogL(i,aLocus,modAlleleFreqs);
  return(logPr);
}

// calculate global logL (all indiv, all loci) with proposed change in allele freqs
double BayesAss::BAgloballogL(vector< vector< vector<double> > > modAlleleFreqs)
{
  double logPr=0.0;
  for(int i=0;i<static_cast<int>(sample.size());i++)
    logPr+=BAlocallogL(i,modAlleleFreqs);
  return(logPr);
}

void BayesAss::PrintAlleleFreqs()
{
  cout << "----------------------------------------------------------------------\n";
  cout << " Posterior Mean Allele Frequency Estimates:                           \n";
  cout << "----------------------------------------------------------------------\n";
  for(int i=0;i<=static_cast<int>(popSet.size())-1;i++)
    {
      cout << "\nPopulation: " << indexTopop[i] << "\n";
      for(int j=0;j<=static_cast<int>(avgAlleleFreqs[i].size())-1;j++)
	{
	  cout << "Locus " << j << ": ";
	  for(int k=1;k<=static_cast<int>(avgAlleleFreqs[i][j].size())-1;k++)
	    cout << intToallele[j][k]<< ":" << avgAlleleFreqs[i][j][k] << " ";
	  cout << "\n";
	}
    } 
  cout << "\n";
}

void BayesAss::BayesAssMCMC(Parameters pvalues, ofstream &outfile, ofstream &logfile)
{
  // initiate random number generator
  gsl_rng *r;	
  r=gsl_rng_alloc(gsl_rng_taus);
  gsl_rng_set (r,pvalues.seed); 

  // acceptance rates
  double afreq_accept_rate=0.0;
  double immanc_accept_rate=0.0;
  double fstat_accept_rate=0.0;
  double migration_accept_rate=0.0;

  long int thinnedSampleCount=1;

  // summary statistics
  vector <vector <vector <double> > > ancestryAssignFreqs;
  for(int i=0; i<static_cast<int>(sample.size()); i++)
    {
      vector <double> myInd;
      for(int j=0; j<=2; j++)
	myInd.push_back(0.0);
      vector < vector<double> > myPop;
      for(int j=0; j<static_cast<int>(popSet.size()); j++)
	myPop.push_back(myInd);
      ancestryAssignFreqs.push_back(myPop);
    }
  vector <vector <double> > fstatSummary;
  vector <double> MeanVar;
  MeanVar.push_back(0.0);
  MeanVar.push_back(0.0);
  for(int i=0; i<static_cast<int>(sample.size()); i++)
    {
      fstatSummary.push_back(MeanVar);
    }
  
  // debugging 
  // PrintAlleleFreqs();

  cout << "\nStarting BayesAss MCMC (Iterations:" << pvalues.MCMCreps << ")\n";
  cout << "----------------------------------------------------------------------\n";
  cout << "Iter:   " << "FrAcc:  " << "ImAcc:  " << "FsAcc:  " << "MgAcc:  " << "logP:      " << "logL:   " << endl;

  logfile << "State \t logL \t logPr";
  for(int i=0;i<static_cast<int>(popSet.size());i++)
    {
      for(int j=0;j<static_cast<int>(popSet.size());j++)
	{
	  logfile << "\t" << indexTopop[i] << "->" << indexTopop[j];
	}
    }

  // time the MCMC run
  Time start;
  
  // initialize local variables for MCMC

  double tempAF[500], tempAP[500];
  double tempAF2[500], tempAP2[500];
  vector< vector< vector<double> > > proposedAlleleFreqs=alleleFreqs; 
  vector< vector< vector<long int> > > proposedMigrantCounts;
  vector< vector< vector<long int> > > zeroMigrantCounts;
  vector<long int> tempRow1;
  vector<vector<long int> > tempRow2; 
  for(int k=0;k<=2;k++)
    tempRow1.push_back(0);
  for(int i=0;i<static_cast<int>(popSet.size());i++)
    tempRow2.push_back(tempRow1);
  for(int i=0;i<static_cast<int>(popSet.size());i++)
    proposedMigrantCounts.push_back(tempRow2);
  zeroMigrantCounts=proposedMigrantCounts;
  
  vector<AugmentIndiv> proposedAugmentedSample;

  double tempMR[500],tempMR2[500];
  double temp2MR[500],temp2MR2[500];
  vector< vector <double> > proposedMigrationRates;

  //begin MCMC

  for(int iter=1;iter<=(pvalues.MCMCreps+pvalues.MCMCburnin);iter++)
    {
      double alpha=0.0;


#ifndef NOALLELEFREQ /* turn off allele freq change proposal for debugging */      
      if(gsl_rng_uniform_pos(r)<=0.2)
	{
	  // propose new allele freqs at a uniformly chosen locus and population
	  long int aLocus=gsl_rng_uniform_int(r,sample[0].genotypes.size());
	  long int aPopulation=gsl_rng_uniform_int(r,popSet.size());
	  proposedAlleleFreqs=alleleFreqs;
	  for(int i=0;i<static_cast<int>(alleleFreqs[0][aLocus].size()-1);i++)
	    {
	      tempAP[i]=alleleFreqs[aPopulation][aLocus][i+1]*pvalues.alleleFreqDirichlet;
	    }
	  gsl_ran_dirichlet(r,alleleFreqs[0][aLocus].size()-1,tempAP,tempAF);
	  
	  // prevent underflows
	  for(int i=0;i<static_cast<int>(alleleFreqs[0][aLocus].size()-1);i++)
	    if(tempAF[i]<10e-8) tempAF[i]=10e-8;
	  
	  for(int i=1;i<static_cast<int>(alleleFreqs[0][aLocus].size());i++)
	    {
	      proposedAlleleFreqs[aPopulation][aLocus][i]=tempAF[i-1];
	    }
	  
	  for(int i=0;i<static_cast<int>(alleleFreqs[0][aLocus].size()-1);i++)
	    {
	      tempAP2[i]=tempAF[i]*pvalues.alleleFreqDirichlet;
	      tempAF2[i]=alleleFreqs[aPopulation][aLocus][i+1];
	    }
      
	  // MH accept/reject
	  alpha=exp(
#ifndef NODATA /* turn off likelihood for debugging */
		    BAgloballogL(aLocus,proposedAlleleFreqs)-
		    BAgloballogL(aLocus)+
#endif /* turn off likelihood for debugging */
		    gsl_ran_dirichlet_lnpdf(alleleFreqs[0][aLocus].size()-1,tempAP2,tempAF2)-
		    gsl_ran_dirichlet_lnpdf(alleleFreqs[0][aLocus].size()-1,tempAP,tempAF));
	  if(gsl_rng_uniform_pos(r)<=alpha)
	    {
	      alleleFreqs[aPopulation][aLocus]=proposedAlleleFreqs[aPopulation][aLocus];
	      afreq_accept_rate=afreq_accept_rate*((iter-1.0)/(double)iter)+1.0/iter;
	    }
	  else
	    afreq_accept_rate=afreq_accept_rate*((iter-1.0)/(double)iter);
	}
#endif /* turn off allele freq change proposal for debugging */

     // propose new immigrant ancestry for a random individual
      proposedAugmentedSample=augmentedSample;
      long int aIndiv=gsl_rng_uniform_int(r,sample.size());
      if(gsl_rng_uniform_pos(r)<=1.0/(2.0*popSet.size()-1.0))
	{
	  proposedAugmentedSample[aIndiv].migrantPopulation=popToindex[sample[aIndiv].sampledPopulation];
	  proposedAugmentedSample[aIndiv].ancestryAge=0;
	}
      else
	{
	  int propPopln=popToindex[sample[aIndiv].sampledPopulation];
	  while(propPopln==popToindex[sample[aIndiv].sampledPopulation])
	    propPopln=gsl_rng_uniform_int(r,popSet.size());
	  proposedAugmentedSample[aIndiv].migrantPopulation=propPopln;
	  if(gsl_rng_uniform_pos(r)<=0.5)
	    proposedAugmentedSample[aIndiv].ancestryAge=1;
	  else
	    proposedAugmentedSample[aIndiv].ancestryAge=2;
	}

      // collect counts from proposedAugmentedSample into proposedMigrantCounts
      proposedMigrantCounts=zeroMigrantCounts;
      for(int i=0;i<static_cast<int>(sample.size());i++)
	++proposedMigrantCounts[popToindex[sample[i].sampledPopulation]][proposedAugmentedSample[i].migrantPopulation][proposedAugmentedSample[i].ancestryAge];
      
      // MH accept/reject
      alpha=exp(
#ifndef NODATA /* turn off likelihood for debugging */
		BAlocallogL(aIndiv,proposedAugmentedSample)-BAlocallogL(aIndiv)+
#endif /* turn off likelihood for debugging */
#ifdef NODATA
		0.0000000001+
#endif
		0.00
#ifndef NOMIGRANTCOUNTPRIOR /* turn off migration count prior for debugging */
		+BAmigratelogP(pvalues,proposedMigrantCounts)-BAmigratelogP(pvalues)
#endif /* turn off migration count prior for debugging */
); 
      if(gsl_rng_uniform_pos(r)<=alpha)
	{
	  augmentedSample=proposedAugmentedSample;
	  for(int i=0;i<static_cast<int>(popSet.size());i++)
	    for(int j=0;j<static_cast<int>(popSet.size());j++)
	      for(int k=0;k<=2;k++)
		migrantCounts[i][j][k]=0;
	  for(int i=0;i<static_cast<int>(sample.size());i++)
	    ++migrantCounts[popToindex[sample[i].sampledPopulation]][augmentedSample[i].migrantPopulation][augmentedSample[i].ancestryAge];
	  immanc_accept_rate=immanc_accept_rate*((iter-1.0)/(double)iter)+1.0/iter;
	}
      else
	immanc_accept_rate=immanc_accept_rate*((iter-1.0)/(double)iter);

#ifndef NOMIGRATIONRATE /* turn off migration rate proposals for debugging */ 
      // propose new migration rate for a randomly chosen column of the migration matrix
      proposedMigrationRates=migrationRates;
      long int aColumn=gsl_rng_uniform_int(r,popSet.size());

      // if Ngens==1 use constraint that sum of off-diagonal rates equals 0.3 
      if(pvalues.Ngens==1)
	{
	  for(int i=0;i<static_cast<int>(popSet.size());i++)
	    {
	      tempMR[i]=migrationRates[i][aColumn]*(pvalues.migPropDirichlet/maxMigration);
	    }
	  gsl_ran_dirichlet(r,popSet.size(),tempMR,tempMR2);
	  
	  // prevent underflows
	  for(int i=0;i<static_cast<int>(popSet.size());i++)
	    if(tempMR2[i]<10e-8) tempMR2[i]=10e-8;
	  
	  for(int i=0;i<static_cast<int>(popSet.size());i++)
	    {
	      proposedMigrationRates[i][aColumn]=maxMigration*tempMR2[i];
	    }
	  
	  for(int i=0;i<static_cast<int>(popSet.size());i++)
	    {
	  temp2MR[i]=tempMR2[i]*pvalues.migPropDirichlet;
	  temp2MR2[i]=migrationRates[i][aColumn]/maxMigration;
	    }
	}

      // else if Ngens==0 use no constraints on rates
      else
	{
	  for(int i=0;i<static_cast<int>(popSet.size());i++)
	    {
	      tempMR[i]=migrationRates[i][aColumn]*(pvalues.migPropDirichlet);
	    }
	  gsl_ran_dirichlet(r,popSet.size(),tempMR,tempMR2);
	  
	  // prevent underflows
	  for(int i=0;i<static_cast<int>(popSet.size());i++)
	    if(tempMR2[i]<10e-8) tempMR2[i]=10e-8;
	  
	  for(int i=0;i<static_cast<int>(popSet.size());i++)
	    {
	      proposedMigrationRates[i][aColumn]=tempMR2[i];
	    }
	  
	  for(int i=0;i<static_cast<int>(popSet.size());i++)
	    {
	  temp2MR[i]=tempMR2[i]*pvalues.migPropDirichlet;
	  temp2MR2[i]=migrationRates[i][aColumn];
	    }
	}

      // MH accept/reject
      alpha=exp(
		#ifndef NODATA
		BAmigratelogP(pvalues,proposedMigrationRates)-BAmigratelogP(pvalues)+
                #endif
		gsl_ran_dirichlet_lnpdf(popSet.size(),temp2MR,temp2MR2)-
		gsl_ran_dirichlet_lnpdf(popSet.size(),tempMR,tempMR2));
      if(gsl_rng_uniform_pos(r)<=alpha)
	{
	  migrationRates=proposedMigrationRates;
	  migration_accept_rate=migration_accept_rate*((iter-1.0)/(double)iter)+1.0/iter;
	}
      else
	migration_accept_rate=migration_accept_rate*((iter-1.0)/(double)iter);
#endif /* turn off migration rate proposals for debugging */ 

#ifndef NOINBREEDING /* turn off inbreeding coefficient proposals for debugging */
      // propose new inbreeding coefficient for a random individual
      proposedAugmentedSample=augmentedSample;
      aIndiv=gsl_rng_uniform_int(r,sample.size());
      proposedAugmentedSample[aIndiv].fstat=
	slidingWindowProposal(gsl_ran_flat(r,-pvalues.windowDelta,pvalues.windowDelta),0,1,augmentedSample[aIndiv].fstat);

      // MH accept/reject
      #ifndef NODATA
      alpha=exp(BAlocallogL(aIndiv,proposedAugmentedSample)-BAlocallogL(aIndiv));
      #endif
      #ifdef NODATA 
      alpha=1.0; 
      #endif
      if(gsl_rng_uniform_pos(r)<=alpha)
	{
	  augmentedSample=proposedAugmentedSample;
	  fstat_accept_rate=fstat_accept_rate*((iter-1.0)/(double)iter)+1.0/iter;
	}
      else
	fstat_accept_rate=fstat_accept_rate*((iter-1.0)/(double)iter);
#endif /* turn off inbreeding coefficient proposals for debugging */


      // print MCMC samples to file BA2.log
      if(((iter % pvalues.MCMCthin)==0)&&(iter>pvalues.MCMCburnin))
	{
	  logfile << "\n" << iter << "\t" << BAmigratelogP(pvalues) << "\t" << BAgloballogL();
	  vector<double> sODR;
	  for(int i=0;i<static_cast<int>(popSet.size());i++)  
	    {
	      double tempSum=0.0;
	      for(int j=0;j<static_cast<int>(popSet.size());j++)
		{
		  if(i!=j) tempSum+=migrationRates[j][i];
		}
	      sODR.push_back(tempSum);
	    }
	  for(int i=0;i<static_cast<int>(popSet.size());i++)
	    {
	      for(int j=0;j<static_cast<int>(popSet.size());j++)
		{
		  if(i==j)
		    {
		      logfile << "\t" << (1.0-sODR[i]);
		    }
		  else
		    {
		      logfile << "\t" << migrationRates[j][i];
		    }
		}
	    }
	}


      // update summary statistics
      if(((iter % pvalues.MCMCthin)==0)&&(iter>pvalues.MCMCburnin))
	{
	  for(int i=0; i<static_cast<int>(sample.size()); i++)
	    for(int j=0; j<static_cast<int>(popSet.size()); j++)
	      for(int k=0; k<=2; k++)
		{
		if((augmentedSample[i].migrantPopulation==j)&&(augmentedSample[i].ancestryAge==k))
		  ancestryAssignFreqs[i][j][k]=ancestryAssignFreqs[i][j][k]*
		    (thinnedSampleCount-1.0)/((double)thinnedSampleCount)+1.0/((double)thinnedSampleCount);
		else
		  ancestryAssignFreqs[i][j][k]=ancestryAssignFreqs[i][j][k]*(thinnedSampleCount-1.0)/((double)thinnedSampleCount);
		}
	  for(int i=0; i<static_cast<int>(sample.size()); i++)
	    {
	      fstatSummary[i][0]=fstatSummary[i][0]*(thinnedSampleCount-1.0)/((double)thinnedSampleCount)+augmentedSample[i].fstat/((double)thinnedSampleCount);
	      fstatSummary[i][1]=fstatSummary[i][1]*(thinnedSampleCount-1.0)/((double)thinnedSampleCount)+(augmentedSample[i].fstat*augmentedSample[i].fstat)/
		((double)thinnedSampleCount);
	    }
	  for(int i=0;i<=static_cast<int>(popSet.size())-1;i++)
	    {
	      for(int j=0;j<=static_cast<int>(avgAlleleFreqs[i].size())-1;j++)
		{
		  for(int k=0;k<=static_cast<int>(avgAlleleFreqs[i][j].size())-1;k++)
		   avgAlleleFreqs[i][j][k]=avgAlleleFreqs[i][j][k]*(thinnedSampleCount-1.0)/((double)thinnedSampleCount)+alleleFreqs[i][j][k]/((double)thinnedSampleCount);
		}
	    }
	  for(int i=0;i<=static_cast<int>(popSet.size())-1;i++)
	    {
	      for(int j=0;j<=static_cast<int>(popSet.size())-1;j++)
		{
		   avgMigrationRates[i][j]=avgMigrationRates[i][j]*(thinnedSampleCount-1.0)/((double)thinnedSampleCount)+migrationRates[i][j]/((double)thinnedSampleCount);
		   m2MigrationRates[i][j]=m2MigrationRates[i][j]*(thinnedSampleCount-1.0)/((double)thinnedSampleCount)+
		     (migrationRates[i][j]*migrationRates[i][j])/((double)thinnedSampleCount);
		}
	    }
	  thinnedSampleCount+=1;
	}

      // print summaries to screen
      if((iter % pvalues.MCMCdisplay)==0)
	{ 
	  cout << showpoint;
	  cout << setw(8) << setfill(' ') << setiosflags(ios::left) << (iter-pvalues.MCMCburnin);
	  cout << setw(6) << setfill('0') << setiosflags(ios::left) << setprecision(3) << afreq_accept_rate;
	  cout << "  ";
	  cout << setw(6) << setprecision(3) << setiosflags(ios::left) << immanc_accept_rate;
	  cout << "  ";
	  cout << setw(6) << setprecision(3) << setiosflags(ios::left) << fstat_accept_rate;
	  cout << "  ";
	  cout << setw(6) << setprecision(3) << setiosflags(ios::left) << migration_accept_rate;
	  cout << "  ";
	  cout << setw(8) << setprecision(6) << setiosflags(ios::left) << BAmigratelogP(pvalues);
	  cout << "   ";
	  cout << setw(10) << setprecision(8) << setiosflags(ios::left) << BAgloballogL(); // << "\n";
	  cout << noshowpoint;
	  if(iter>pvalues.MCMCburnin)
	    {
	      Time current;
	      cout << "  " << floor(((current.delta(&start)/(double)iter)*(pvalues.MCMCreps+pvalues.MCMCburnin-iter))/60) << ":";
	      cout << setw(2) << setprecision(2) << setfill('0') << (int)(((current.delta(&start)/(double)iter)*(pvalues.MCMCreps+pvalues.MCMCburnin-iter))) % 60 << "\n";
	    }
	  else
	    {
	      cout << endl;
	    }
	}	  
      // print individual population assignments to output file
      outfile << "\n iter: " << iter;	  
      for(int i=0;i<static_cast<int>(sample.size());i++)
	{
	  outfile << setw(4) << augmentedSample[i].migrantPopulation << ":" << augmentedSample[i].ancestryAge;
	}
    }
  gsl_rng_free(r);

  // print population assignment freqs and fstats to screen
  for(int i=0;i<static_cast<int>(sample.size());i++)
    {
      cout << "\n Individual: " << indexTolabel[i] << " Population: " << sample[i].sampledPopulation << "\n";
      for(int j=0; j<=2; j++)
	{
	  cout << " Age=" << j << "  ";
	  for(int k=0; k<static_cast<int>(popSet.size()); k++)
	    {
	      cout << indexTopop[k]+":" << setw(6) << setprecision(3) << showpoint << ancestryAssignFreqs[i][k][j] << "  ";
	      cout << noshowpoint;
	    }
	  cout << "\n";
	}
      cout << " Fstat: "  << setw(5) << setprecision(3) << fstatSummary[i][0];
      if((fstatSummary[i][0]-1.96*sqrt(fstatSummary[i][1]-fstatSummary[i][0]*fstatSummary[i][0]))<0)
	cout << " (" << showpoint << setw(5) << setprecision(3) << 0.0000 << "," << noshowpoint;
      else
	cout << " (" << setw(5) << setprecision(3) << fstatSummary[i][0]-1.96*sqrt(fstatSummary[i][1]-fstatSummary[i][0]*fstatSummary[i][0]) << ",";
      if((fstatSummary[i][0]+1.96*sqrt(fstatSummary[i][1]-fstatSummary[i][0]*fstatSummary[i][0]))>1.0)
	cout << setw(5) << showpoint << setprecision(3) << 1.0000 << ")\n" << noshowpoint;
      else
	cout << setw(5) << setprecision(3) << fstatSummary[i][0]+1.96*sqrt(fstatSummary[i][1]-fstatSummary[i][0]*fstatSummary[i][0]) << ")\n";
    }

  // print mean posterior allele freqs
  PrintAlleleFreqs();


  /*  vector<double> sumODR;
  for(int i=0;i<static_cast<int>(popSet.size());i++)  
    {
      double tempSum=0.0;
      for(int j=0;j<static_cast<int>(popSet.size());j++)
	{
	  if(i!=j) tempSum+=avgMigrationRates[j][i];
	}
      sumODR.push_back(tempSum);
      } */

  // print migration rates to screen
  cout << "\n----------------------------------------------------------------------\n";
  cout << "Posterior Distribution of Migrant Proportions Among Populations\n";
  cout << "----------------------------------------------------------------------" << endl;
  cout << "       E(mij) [untransformed]  E(mij)=1 [normalized]";
  cout << setfill(' ');
  for(int i=0;i<static_cast<int>(popSet.size());i++)
    {
       for(int j=0;j<static_cast<int>(popSet.size());j++)
	{
	      cout << "\n" << " " << indexTopop[j] << "->";
	      cout << indexTopop[i] << ": " << setw(8) << setprecision(5) << setiosflags(ios::left) << avgMigrationRates[j][i];
	      if((avgMigrationRates[j][i]-1.96*sqrt(m2MigrationRates[j][i]-avgMigrationRates[j][i]*avgMigrationRates[j][i]))<0)
		cout << " (" << showpoint << setw(5) << setprecision(3) << 0.0000 << "," << noshowpoint;
	      else
		cout << " (" << setw(5) << setprecision(3) << avgMigrationRates[j][i]-1.96*sqrt(m2MigrationRates[j][i]-avgMigrationRates[j][i]*avgMigrationRates[j][i])  << ",";
	      if((avgMigrationRates[j][i]+1.96*sqrt(m2MigrationRates[j][i]-avgMigrationRates[j][i]*avgMigrationRates[j][i]))>1.0)
		cout << setw(5) << showpoint << setprecision(3) << 1.0000 << ")" << noshowpoint;
	      else
		cout << setw(5) << setprecision(3) << avgMigrationRates[j][i]+1.96*sqrt(m2MigrationRates[j][i]-avgMigrationRates[j][i]*avgMigrationRates[j][i]) << ")";

	      cout << "  " <<  setw(8) << setprecision(5) << setiosflags(ios::left) << avgMigrationRates[j][i]*popSet.size(); // normalized
	      if((avgMigrationRates[j][i]-1.96*sqrt(m2MigrationRates[j][i]-avgMigrationRates[j][i]*avgMigrationRates[j][i]))<0)
		cout << " (" << showpoint << setw(5) << setprecision(3) << 0.0000 << "," << noshowpoint;
	      else
		cout << " (" << setw(5) << setprecision(3) << popSet.size()*avgMigrationRates[j][i]-1.96*sqrt(m2MigrationRates[j][i]-avgMigrationRates[j][i]*avgMigrationRates[j][i])*popSet.size()  << ",";
	      cout << setw(5) << setprecision(3) << popSet.size()*avgMigrationRates[j][i]+1.96*sqrt(m2MigrationRates[j][i]-avgMigrationRates[j][i]*avgMigrationRates[j][i])*popSet.size() << ")";
	}
    }
  
  cout << "\n";
  cout << "----------------------------------------------------------------------\n";
  
  Time current;
  cout << "Elapsed Time: " << current.delta(&start) << " seconds\n"; 

}
